// ml:run = time -p $bin
#include <iostream>
#include <fstream>
#include <iomanip>
#include <numeric>
#include <vector>
#include <cmath>
#include <Eigen/SparseCholesky>
#include <Eigen/SparseLU>
#include <boost/math/constants/constants.hpp>

namespace hexagon
{
    using SpMat = Eigen::SparseMatrix<double>;
    using triplet = Eigen::Triplet<double>;

    double constexpr pi = boost::math::constants::pi<double>();
    double constexpr eps = 1e-6;

    // number of all points
    int div{1};
    int n{15};
    int width{5};
    int height{3};
    double length{2. / std::sqrt(3)};
    double delta{length / 2. / div};
    double g{-5. / 4.};

    // f values, counterclockwise from degree 0
    std::vector<double> f{0, 1, 2, 3, 4, 5};
    // std::vector<double> f{4, 4, 4, 4, 4, 4};
    std::vector<double> cc{0, 0, 0, 0, 0, 0};

    struct point { double x; double y; };

    template <class T>
    T abs(T x) { return x < 0 ? -x : x; }

    // judge if x in [l, r]
    bool inrange(int x, int l, int r) { return l <= x && x <= r; }

    bool boundary(int x, int y)
    {
        if (y == 0          && inrange(x, div, 3 * div)) return true;
        if (y == height - 1 && inrange(x, div, 3 * div)) return true;
        if (div - x     == y || div - x     == height - 1 - y) return true;
        if (x - 3 * div == y || x - 3 * div == height - 1 - y) return true;
        return false;
    }

    bool inside(int x, int y)
    {
        if (y <= div) return inrange(x, div - y + 1, 3 * div + y - 1);
        auto d = height - 1 - y;
        return inrange(x, div - d + 1, 3 * div + d - 1);
    }

    bool outside(int x, int y)
    {
        return !boundary(x, y) && !inside(x, y);
    }

    point real_pos(int x, int y)
    {
        int cx{2 * div};
        int cy{div};
        int dx{x - cx};
        int dy{y - cy};
        return {dx * delta, dy * std::sqrt(3) * delta};
    }

    auto identity(int x, int y) { return y * width + x; }

    double value(int x, int y)
    {
        auto p = real_pos(x, y);
        if (p.x == 0 && p.y == 0) return f[0];
        double cos{p.x / std::hypot(p.x, p.y)};
        double angle{std::acos(cos)};
        if (p.y < 0) angle = 2. * pi - angle;
        int id = angle / (pi / 3.) + eps;
        cc[id]++;
        return f[id];
    }

    namespace impl
    {
        // width  = 4 * div + 1
        // height = 2 * div + 1
        void init(SpMat & A, Eigen::VectorXd & b, std::vector<triplet> & coe)
        {
            // auto id = identity(3 * div, 2 * div);
            // coe.emplace_back(id, id, 1);
            // b[id] = 0;
            // std::cerr << "set zero point: " << id << "\n";

            int out_count{0};
            int in_count{0};
            int on_count{1};

            for (auto j = 0; j < height; j++)
            for (auto i = 0; i < width; i++) {
                // if (i == 3 * div && j == 2 * div) continue;

                auto id    = identity(i, j);
                auto up    = identity(i, j + 1);
                auto down  = identity(i, j - 1);
                auto left  = identity(i - 1, j);
                auto right = identity(i + 1, j);

                if (boundary(i, j)) {
                    on_count++;
                    // top
                    if (j == height - 1) {
                        coe.emplace_back(id, id,    +8);
                        coe.emplace_back(id, right, -3);
                        coe.emplace_back(id, left,  -3);
                        coe.emplace_back(id, down,  -2);
                        b[id] = 3 * value(i, j) * delta * delta + 2 * std::sqrt(3) * g * delta;
                        continue;
                    }
                    // bottom
                    if (j == 0) {
                        coe.emplace_back(id, id,    +8);
                        coe.emplace_back(id, right, -3);
                        coe.emplace_back(id, left,  -3);
                        coe.emplace_back(id, up,  -2);
                        b[id] = 3 * value(i, j) * delta * delta + 2 * std::sqrt(3) * g * delta;
                        continue;
                    }

                    auto dx = i > 2 * div ? +1 : -1;
                    auto dy = j >     div ? +1 : -1;
                    auto tid1 = identity(i, j - dy);
                    auto tid2 = identity(i - dx, j);
                    coe.emplace_back(id, id,   +8);
                    coe.emplace_back(id, tid1, -2);
                    coe.emplace_back(id, tid2, -6);
                    b[id] = 3 * value(i, j) * delta * delta + 4 * std::sqrt(3) * g * delta;
                    continue;
                }
                if (outside(i, j)) {
                    out_count++;
                    coe.emplace_back(id, id, 1);
                    b[id] = 0;
                    continue;
                }

                // inside
                in_count++;
                coe.emplace_back(id, id,    +8);
                coe.emplace_back(id, right, -3);
                coe.emplace_back(id, left,  -3);
                coe.emplace_back(id, up,    -1);
                coe.emplace_back(id, down,  -1);

                b[id] = 3 * value(i, j) * delta * delta;
            }

            std::cerr << "on_count = " << on_count << "\n";
            std::cerr << "in_count = " << in_count << "\n";
            std::cerr << "out_count = " << out_count << "\n";
            std::cerr << "tot count = " << on_count + in_count + out_count << "\n";

            for (auto i : cc) std::cerr << i << " ";
            std::cerr << "\n";

            // std::cerr << "size = " << coe.size() << "\n";
            // for (auto i : coe)
            //     std::cerr << i.row() << " " << i.col() << " " << i.value() << "\n";

            A.setFromTriplets(coe.begin(), coe.end());
        }
    }

    // hexagon edge length is 2/sqrt(3)
    void solver(int division = -1)
    {
        // division for half an edge
        if (division == -1) {
            std::cout << "please enter a division number of a half length edge.\n";
            std::cin >> division;
        }

        div = division;
        delta = (length / 2.) / div;
        g = -std::accumulate(f.begin(), f.end(), 0.) / 12.;
        width  = 4 * div + 1;
        height = 2 * div + 1;

        // total points in a rectangular area for convenient. Actually,
        // 6 * div * div + 4 * div + 1 total points
        n = (2 * div + 1) * (4 * div + 1);
        auto actual = 6 * div * div + 4 * div + 1;

        std::cerr << "div = " << div << "\n";
        std::cerr << "n = " << n << "\n";
        std::cerr << "actual = " << actual << "\n";
        std::cerr << "g = " << g << "\n";
        std::cerr << "width = " << width << "\n";
        std::cerr << "height = " << height << "\n";

        std::vector<triplet> coefficients;
        coefficients.reserve(5 * n);

        SpMat A(n, n);
        // A.reserve(n * n + n);
        Eigen::VectorXd b(n);

        // initialize A from bottom to top
        impl::init(A, b, coefficients);

        // std::fstream ffout{"data.in"};
        // ffout << n << "\n";
        // for (int i = 0; i < n; i++) {
        //     for (int j = 0; j < n; j++) {
        //         ffout << A.coeff(i, j) << " ";
        //     }
        //     ffout << b[i] << "\n";
        // }


        Eigen::VectorXd x(n);

        // std::fstream fin{"data.out"};
        // for (int i{0}; i < n; i++)
        //     fin >> x[i];

        // Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
        // solver.compute(A);

        Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
        A.makeCompressed();
        solver.analyzePattern(A);
        solver.factorize(A);
        if (solver.info() != Eigen::Success) {
            std::cerr << "decomposition failed.\n";
            return;
        }
        x = solver.solve(b);
        if (solver.info() != Eigen::Success) {
            std::cerr << "solving failed.\n";
            return;
        }

        // std::cerr << "solution size: " << x.size() << "\n";
        // for (int i{0}; i < n; i++)
        //     // fin << x[i] << "\n";
        //     std::cerr << x[i] << "\n";

        std::ofstream fout{"points.dat"};
        // output to tecplot data format
        fout << "Title=\"solution\"\nVariables=\"x\", \"y\", \"u\"\n";
        for (int px = 0; px < width; px++)
            for (int py = 0; py < height; py++) {
                if (outside(px, py)) continue;
                auto p = real_pos(px, py);
                fout << std::fixed << std::setprecision(6)
                    << p.x << " " << p.y << " " << x[identity(px, py)] << "\n";
            }


        double eps = 0;
        std::vector<double> sum(n);
        for (auto i : coefficients) {
            int row = i.row();
            int col = i.col();
            double v = i.value();
            sum[row] += x[col] * v;
        }
        for (int i{0}; i < n; i++)
            eps = std::max(eps, abs(sum[i] - b[i]));
        std::cerr << "max eps is: " << eps << "\n";
    }
}

int main()
{
    hexagon::solver(40);
}

